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Resolving infall caustics in dark matter halos 
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We have found that the phase-space of a dark matter particles assembling a galactic halo in cosmological 
N-body simulations has well defined fine grained structure. Recently accreted particles form distinctive ve- 
locity streams with high density contrast. For fixed observer position these streams lead to peaks in velocity 
i— i distribution. Overall structure is close to that emerging in the secondary infall model. 
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1. Introduction. 

o 

Precise knowledge of phase-space distribution of the Galactic halo dark matter particles is important for the 
<X3 dark matter search experiments. Phase-space model, which is idealized but in certain aspects more realistic 

1 ^ i than e.g. traditional "isothermal" model and capable to describe fine-grained structure of the phase-space, was 

constructed in Refs. [1, 2]. A distinctive feature of corresponding particle distributions is that the highest energy 
T ^ particles have discrete values of velocity [3]. Qualitatively this conclusion is based on the following observation. 

Initial velocity dispersion of cold dark matter particles is very small, therefore one can consider that they occupy 
f^*) thin 3-dimensional sheet in a 6-dimensional phase space, v = Hf. Later on, when structure forms, the Hubble 

law would remain intact at large distances for isolated halos, but inside halos this sheet rolls and wraps around 

itself forming growing overdensity. Owing to the Liouville theorem the occupied phase-space sheet can neither 
<0> cross itself nor be perforated. Resulting phase-space of the halo dark matter particles is depicted in Fig. 1. 

£SJ The model can be solved exactly with the assumption of self-similar, radial infall [4, 5]. In Ref. [2] the model 
was extended to include angular momentum in a self-similar way. (Phase-space shown in Fig. 1 corresponds 
to zero angular momentum.) With angular momentum included, the model explains not only flat rotational 
curves at intermediate distances, but also p cx r~ Q behavior of density, with a close to 1 in the cusps [2]. Note 
■ that assumptions of self-similarity and spherical symmetry are needed only to get exact solutions. With these 
C3 assumptions relaxed, main qualitative features of the model would survive. In phase-space we would still have 
dense separated folds occupied by matter. Without self-similarity fold separations will change, while without 
spherical symmetry folds will became triaxial, but for a number of phenomenological applications this is inessen- 
tial. For example, axion detectors have very good energy resolution and if dark matter is made of axions, such 
fine structure of velocity space can be resolved [2] . 

Near the turn-around radius for each fold the density of dark matter in configuration space increases, such 
regions of space are called "caustics". It has been suggested that caustics would increase the dark matter anni- 
hilation rate [6, 7, 8] and that they can be probed by properly stacking the weak-lensing signal of a reasonable 
number of dark halos [9]. An evidence for inner caustic rings distributed according to the predictions of the 
self-similar infall model has been found in the Milky Way [10, 11, 12], in other isolated spiral galaxies [10, 13], 
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Figure 1. Phase space distribution of the halo dark matter particles at a fixed moment of time in the infall 
model. The solid curves represent occupied phase-space cells. The vertical dashed line corresponds to the ob- 
server position. Each intersection of the solid and dashed lines corresponds to a velocity peak in the velocity 
distribution measured by an observer. A small scale sub-clump falling into the galaxy for the first time is also 
shown. Adapted from Ref. [2]. 



and in galaxy clusters [14], while there is clear evidence of continuing infall in our Local Group [15, 16] and also 
in the ensemble of halos in a wide range of scales from small galaxies to galaxy clusters [17]. 

Non zero primordial (e.g. thermal) velocities of dark matter particles can be accounted for analytically [7]. 
Their effect is negligibly small as far as dark matter is cold. In Ref. [18] a formalism was developed which allows an 
analytical treatment of phase space streams and caustics without assumptions of spherical symmetry, continuous 
or smooth accretion, and self-similar infall for the formation of dark matter halos. However, the overall problem 
is difficult. Such infall picture will be destroyed e.g. by a major recent merger event. Influence of steady infall of 
small-scale clumps, as the one shown in Fig. 1 and which is outside of frameworks of self-similar infall model, is 
not clear also. Regarding influence of sub-clumps we would like to make several comments: i) Major fraction of 
dark matter resides outside of clumps and volume filling factor occupied by clumps is small. Outside of clumps 
the phase space may correspond to the infall model, ii) Clumps themselves follow infall trajectories in phase 
space. Inside host halo the clumps are tidally disrupted eventually. Released dark matter particles will continue 
to follow infall trajectories. The velocity dispersion in the resulting tidal streams may be smaller than in the 
clumps themselves since particles are released near the clump boundary were relative velocity is small, iii) The 
number of the observed dwarfs in our Galaxy halo is orders of magnitude less then the numerical simulations of 
ACDM model predict. A solution to this problem can be found on the assumption that the primordial power on 
small scales is reduced. Realization of this scenario in nature would bring realistic situation closer to an idealized 
infall model. 

Ultimately, one would like to resolve fine details of phase-space structure of collapsing halos in direct N-body 
approach. Already in early simulations [19] it was observed that the collapse of galactic halos in the cosmolog- 
ical setup proceeds gently and not via violent relaxation. A strong correlation of the final energy of individual 
particles with the initial one has been found indicating overall validity of infall picture. However, detection of 
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caustic surfaces in N-body numerical modeling is a challenging problem. To understand and resolve the problem 
of fine structure of dark matter phase-space numerically, a dedicated very high resolution numerical simulations 
are needed. (To tackle the issue, novel numerical techniques are also emerging, superseding simple increase of the 
number of particles [20, 21, 22].) Indeed, in simulations each "particle" exceeds by many orders of magnitude the 
solar mass, which is 60 to 80 orders of magnitude larger than the mass of dark matter candidates being simulated. 
This leads to artificial two body scattering, energy transfer and to distorted orbits, see e.g. [23], while in reality 
the dark matter particles are collisionlcss and pass unperturbed past each other. Early simulations contained 
only a few thousand particles, and realistic study of phase space was hopeless. Kincmatically cold infall streams 
were observed in more recent high resolution simulations. In Ref. [24] it was found that the motions of the most 
energetic particles were strongly clumped and highly anisotropic. However, these simulations still produce rather 
smooth mass profiles without high density caustics and the mean contribution of an individual stream to local 
dark-matter distribution near the Sun position in the Galaxy was found to be of order f n 0.3% . This value of 
f n for the solar neighborhood is actually in agreement with predictions of self-similar infall model for e = 1, see 
Ref. [2], where parameter e defines the shape of initial overdensity, 5M(ri)/M(ri) — (Mo / M (ri)) e . 

Until recently, phase space picture qualitatively similar to the one presented in Fig. 1 was not reported from 
N-body simulations, see e.g. Ref. [25]. However, in Ref. [26] it was shown that particle orbits in simulated cosmo- 
logical cold dark matter halos are surprisingly regular and periodic, while the phase-space structure of the outer 
halo shares, qualitatively, some of the properties of the classical self-similar secondary infall model. Namely, the 
positions of six outermost resolved caustic were in very good agreement with the infall model for e = 1. The ratios 
of caustic radii to their turnaround radii also agreed quite well with the model. However, reported trajectories 
in the radial velocity - radius plane were severely broadened producing only small (< 10%) enhancements in the 
spherical density profile. 

In this paper we aim to reanalyze the phase-space structure of the outer halo. Our main point is based on 
the following observation. Realistic halos are triaxial and mass accretion is non-radial. In this situation the 
averaging over all directions in configuration space (as usually done) will not produce meaningful representation 
of the phase-space structure of two dimensional (v r ,r) subspace. Mapping in this way, say, even infinitely thin 
ellipsoidal surfaces, will produce broad band distribution and real nature of phase-space sheets may be lost. To 
reveal it, one should try to display phase space with as small range for averaging, as statistics allows. To verify 
this proposition we construct phase-space by doing averages in a sequence of decreasing solid angles around some 
arbitrarily chosen direction in configuration space. If phase-space becomes sharper in this sequence, then our 
proposition is correct. 

2. Numerical sumulation. 

We analyzed a N-body simulation of a galaxy size halo for which the initial conditions were generated by 
Stoehr et al., Ref. [27], adopting a LCDM cosmology, with fi ro = 0.3, tt L = 0.7, H = 70 km s" 1 Mpc" 1 and 
og = 0.9. This "Milky Way" like halo was selected from an intermediate-resolution simulation as a halo with 
a maximum rotational velocity at z = approximately equal to the Milky Way value as a relatively isolated 
halo, which suffered its last major merger at z > 2. We used the highest resolution initial conditions of that 
set, which has a dark matter particle mass of mdm — 2 x 10 5 M //i, so that the halo at z — was resolved by 
more than 11 million particles within the viral radius. The total mass of the halo withing this radius therefore 
is M v i r — 2.3 x 10 12 Mq/H. The simulation was performd using the cosmological N-Body code Gadget [28, 29] 
with the gravitational softening set to 0.2 kpc/h. 

3. Discussion. 

Phase-space of dark matter particles is 6-dimentional, while we have only 2 dimensions to display the results. 
We start with presenting it in (v r ,r) subspace. In general non-spherically symmetric situation the occupation 
pattern of this subspace will depend upon direction chosen from the halo center. Having limited number of 
particles we cannot fix direction to a particular value, however. To accumulate statistics we simply integrate over 
all directions, as it is usually done. The result is shown in Fig. 2. We see host halo extending to about 300 kpc 
in radius and to 400 km/s in velocity. We also see a large number of infalling sub-clumps, but no strong evidence 
for infall trajectories. Even particles falling in for the first time are forming very wide band in the portion of 
phase-space between 400 kpc and 1 Mpc. Reason for this is clear: stacking even infinitely thin trajectories of the 
sort presented in Fig 1 will result in their dilution in overall ensemble if trajectories differ from each other. 
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Figure 2. Result of N-body simulation. Projection of phase-space of a halo dark matter particles to (r, v r ) plane, 
being averaged over all angles in configuration space. 
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Figure 3. The same as in Fig. 2 but now averaging is restricted to solid angle with opening 6 = 40°. Dotted line 
represents "an observer" position at 150 kpc where the velocity distributions of Fig. 4 are constructed. 



To alleviate this problem we can sum up trajectories only in a limited range of viewing directions. To realize 
this procedure we choose arbitrary reference direction and sum occupation numbers of phase-space cells only if 
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Figure 4. Velocity distribution function at the observer position r = 150 kpc. Averages are done for fixed 
Ar = 2 kpc and solid angles with openings 8 = 40° (thick solid line) and 9 = 20° (thin dotted line). 

direction to a given cell lies within the cone with opening angle 9 around the reference direction. The result is 
shown in Fig. 3 for 9 — 40°. 1 ' Phase-space becomes much sharper now. First infall trajectory can be identified 
and several folds similar to those in Fig. 1 can be visually recognized in the host halo. This tendency continues 
with decreasing solid angle. 

To quantify the tendency we constructed velocity distributions at fixed observer position at r = 150 kpc for 
different values of 9, these are presented in Fig. 4. Thick solid line represents averaging inside solid angle with 
opening 9 — 40° while thin dotted line corresponds to 9 = 20°. Peaks in the velocity distributions correspond 
to crossings of velocity streams in phase-space, see Fig. 3. While amplitudes of peaks are subject to Poisson 
fluctuations to some extent, and occasionally an amplitude may be influenced by the presence of a sub-clump in 
the volume of averaging, overall, these peaks and their average amplitude are real. This can be verified e.g. by 
comparing positions of peaks in Fig. 4 with streams in Fig. 3 and by changing observer position, which we have 
also done. 

Notably, first and second peaks, at v ~ —400 km/s and v w 400 km/s correspondingly, are well separated 
and isolated from the background already at 9 = 40° . We also see that with decreasing 9 the average amplitude 
of peaks becomes consistently higher and they become more narrow. On the average, the peak amplitudes rise 
by a factor of few. This can be compared with small 10% density contrast which was found in Ref. [26] and 
which we also observe at 9 = Air when averaging is done over all angles. This proves our proposition - "blurring" 
of phase-space at large r in the current N-body simulations is due to averaging, which is extraneous to many 
phenomenological applications. In reality the phase space is sharper but we cannot evaluate at present to what 
extent. When we diminish the solid angle of integration even further, peaks continue to grow, but statistics 



1 )" Hairy" tail in the region of distances form 300 kpc to 700 kpc may be a numerical artifact. Similar distortions appear in idealized 
spherically symmetrical situation without initial small-scale perturbations when the problem is solved using 1-dimentional N-body 
simulations instead of methods developed in Refs [4, 5]. If true, this indicates that the processes of two body scattering are distorting 
phase-space in our simulation. 
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Figure 5. The same as in Fig. 3. However, now |?7| ■ sign(v r ) insted of v r is used. 

decreases and noise starts to appear. We cannot reach the regime when the amplitude of the peaks becomes 
^-independent, this will require more extensive numerical simulation. 

However, the phase-space becomes significantly sharper already in the current simulation if instead of v r the 
velocity modulus is chosen to display a two dimensional projection of the six dimensional phase space. In Fig. 5 
we present this by plotting |u| multiplied by the sign of v r . Multiplication by sign of v r is needed to separate 
incoming and outcoming velocity streams. In this figure several largest velocity streams are separated from the 
"noise" completely. They stay separated and well defined up to a small distances of order of the Sun position in 
the Galaxy. Relative mass fraction in these streams is small however. The phase-space is blurred for the earlier 
folds. Again, to clarify the situation in these region one needs more extensive numerical simulation. 

Why the phase-space is sharper in \v\ as compared to v r l The reason can be that v r is simply not a suitable 
coordinate in generic non-sphcrical situation. But whatever the reason, for a number of applications it is enough 
to show that a choice of coordinates exists where phase streams are well defined and separated even after averag- 
ing. This is relevant, for instance, for direct dark matter detection experiments, in particular for axion searches, 
where |u| is important, but not v r . 

We would like to comment also on the following peculiarity of Fig. 5 related to the empty band which is 
passing through the middle of the plot along line \v\ = 0. It appears because \v\ cannot be smaller than the value 
of the transverse component of the velocity near the turnaround point. Therefore, this portion of the phase-space 
is unoccupied, contrary to v r projection. The boundary of this band immediately tells us the value of transverse 
velocity and consequently the value of angular momentum j. The fact that boundary of the band stays parallel 
to \v\ = axis tells that j(r) oc r. 
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4. Conclusions. 

We have demonstrated that the phase-space of a galactic halo emerging in cosmological N-body simulations 
has sharp fine grained structure. The density contrast is high at velocity streams, at least an order of magnitude 
higher than previously reported. This is true for at least several latest streams, i.e. for those which correspond 
to recently accreted particles. Overall structure is close to the one emerging in the infall model. Our results can 
be important for the dark matter search experiments. 
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